(* SIMULATION OF THE JLV : Mirlees with extensive margin         *)
(* Version from May, 2009  Benthamite S.W.F.        8 JULY 2010 *)
(* (c) by Jacquet - Lehmann - Van der Linden (JLV) ou (BEL)      *)
ClearAll["Global`*"];
Needs["PlotLegends`"];

(* Calibrated Parameters *)
ϵ=0.25;                         (* Elasticity of labor supply along a linear tax schedule *)
σ=0.8;                               (* Concavity of individuals' preferences *)
jmax=2000;                         (* number of inside loops 200/16 *)     
Jmax=1500;                         (* Jmax 100/12 *)
kmax=25;                          (* number of medium loops 25/2 *)      
lmax=20;                         (* number of outside loops 20/2 *)    
paralissage=2;   (* smoothing parameter. The higher, the smoother the estimated distribution, but the more fare away from the observations. Benchmark value is 1 *)
TL=1;   (* Length of the time period year=1, month = (1/12), week = (1/52), day=(1/365.25) *)
Pareto=1/2.2;         (* Value of the parameter ((1-F)/(Y*f)) at the upper end of the EARNINGS distribution in the actual economy  Value for 2012:  *)Paretoskill=(Pareto/(1+ϵ)); (* Value of the parameter ((1-F)/(Y*f)) at the upper end of the SKILL distribution in the actual economy  *) 
Employ[w_]:=0.7+0.1*((((w-wmin)/(wmax-wmin)))^0.3);     (*  Approximation of the skill-specific employment rate in the current economy *)
Elastbottom=0.5;
Elasttop=0.4;
Elast[j_]:=Elasttop+(Elastbottom-Elasttop)*(((jmax-j)/jmax))^5;     (*  Approximation of the elasticity of the employment rate in the current economy wrt to ... *)
MarginTaxActualfix=0.108; (*0.279 *)
TaxActual[Y_]:=-841.4+MarginTaxActualfix*Y;(* -4024.9*TL+MarginTaxActualfix*Y;  -841.4*TL+ *)  (*  Approximation of the tax function in the current economy, for singles with no kids *)
bactual=12*(475+237)*TL;(*2381; *)(* 475€: RSA mensuel en 2012, 237€: montant moyen des APL *)
MarginTaxActual[Y_]:=MarginTaxActualfix;       (*  Approximation of the marginal tax rate in the current economy *)
(* Intensive utility function *)
(* The V function *)
V[U_]:=If[σ==1,Log[U+1] ,(((U+1)^(1-σ))/(1-σ))];(* The +1 is here to keep finite marginal utility of consumption when b=0  *)
Vp[U_]:=(U+1)^(-σ); 
Vs[U_]:=(-σ)*((U+1)^(-σ-1));
InvV[x_]:=If[σ==1,Exp[x],((1-σ)x)^((1/(1-σ)))]-1;   
(* The W function *)
W[C_,Y_,w_]:=C-((((Y/w)))^(1+(1/ϵ)));       (* utility levels *)
Γ[U_,Y_,w_]:=InvV[U]+((((Y/w)))^(1+(1/ϵ)));  (* consumption as a function of utility levels *)
Ww[C_,Y_,w_]:=(1+(1/ϵ))*((Y^(1+(1/ϵ))))*((w^(-2-(1/ϵ))));  (*  marginal utility of skill *)
WC[C_,Y_,w_]:=1;     (*  marginal utility of consumption *)
WCw[C_,Y_,w_]:=0;     (*  cross-second derivative of consumption and skill *)
WY[C_,Y_,w_]:=-(1+(1/ϵ))*((Y^(1/ϵ)))*((w^(-1-(1/ϵ))));     (*  marginal utility of earnings *)
U[C_,Y_,w_]:=V[W[C,Y,w]];
UY[C_,Y_,w_]:=Vp[W[C,Y,w]]*WY[C,Y,w];
UC[C_,Y_,w_]:=Vp[W[C,Y,w]]*WC[C,Y,w];
Uw[C_,Y_,w_]:=Vp[W[C,Y,w]]*Ww[C,Y,w];
UCw[C_,Y_,w_]:=Vp[W[C,Y,w]]*WCw[C,Y,w]+Vs[W[C,Y,w]]*Ww[C,Y,w]*WC[C,Y,w];
YLaissezFaire[w_]:=(w^(1+ϵ)) * (((((1+(1/ϵ))))^(-ϵ)));     (*  Earnings in the absence of taxes *)
Skill[Y_,T_,Tm_]:=((Y^(1/(1+ϵ))))*((((((1+(1/ϵ))/(1-Tm))))^(((ϵ/(1+ϵ))))));     (* Inverse the f.o.c of the intensive program to give skill level as a function of earnings Y, level of Tax T and marginal tax rate Tm *)

data=Import["revdisp_i.xls",Path->{"\\\\Vboxsvr\\dd\\Économie\\Données\\ERFS_2012"}]; (*weekly_revdisp_i *)
Nmax=1000;(* Number of observations in the original data 2353 / 760  / 598*)
Nmin=1;  (* 1 Number of deleted observations corresponding to too low earnings 8  1*)

Kernel[x_]:=If[-1<x≤1,0.75(1-(x^2)),0];   (* Kernel function *)

(* data extraction *)

For[j=Nmin,j≤Nmax,j=j+1,
 Yactual[j]=Extract[data,{1,j,1}]*TL;    (* 1.1211* For version 7 of Mathematica 1000-943.31+  *52 *)
skill[j]=Skill[Yactual[j],TaxActual[Yactual[j]],MarginTaxActual[Yactual[j]]];
];

HH[w_,taille_]:=Sum[Kernel[((w-skill[j])/taille)],{j,Nmin,Nmax}];     (*  Gives a Kernel estimation of the disitribution of earnings among employed, individuals, i.e. h=K*f, in the current economy  *)
MMean=((Sum[skill[j],{j,Nmin,Nmax}])/(Nmax-Nmin+1));
MMean2=((Sum[(((skill[j]))^2),{j,Nmin,Nmax}])/(Nmax-Nmin+1));  
taille=paralissage*((((MMean2-((((MMean))^2))))^(1/2)))*(2.34/(((Nmax-Nmin+1))^(1/5))) ;
wmax=skill[Nmax]*(jmax/Jmax);                    (* highest skill level *)
wmin=(wmax/(jmax*1000))(*skill[Nmin]*);                                   (* lowest skill level *)
dw=((wmax-wmin)/(jmax+1));
(*  Participation decisions  *)

P[u_,b_,w_]:=u-U[b,0,w];                (* Threshold function sorting participating against nonemployed individuals of the same skill level *)

(*  calibration from an approximation of the current economy *)

For[j=0,j≤jmax,j=j+1,
w[j]=wmin+((j*(wmax-wmin))/jmax);
ff[j]=HH[w[j],taille];   (* unormalized approximation of the unconditionnal density of skill w, estimated by a Kernel procedure, before the extension by Pareto distribution for the upper tail, indexed by j  *)
YLF[j]=YLaissezFaire[w[j]];
Ycurrent[j]=w[j]*(((((ϵ/(1+ϵ))*w[j]*(1-MarginTaxActualfix)))^ϵ));
(* gross earnings in the current economy at skill w[j] *)
Ucurrent[j]=U[Ycurrent[j]-TaxActual[Ycurrent[j]],Ycurrent[j],w[j]];

β [j]=Elast[j]/((Ycurrent[j]-TaxActual[Ycurrent[j]]-bactual)*UC[Ycurrent[j]-TaxActual[Ycurrent[j]],Ycurrent[j],w[j]]*(1-Employ[w[j]]));                                                                                   (* Parameter for the conditionnal distribution of χ adjusted to match a reasonnable elasticity of participation in the current economy *)
α[j]=β[j]*P[Ucurrent[j],bactual,w[j]]+Log[((1-Employ[w[j]])/(Employ[w[j]]))];                                                     (* Parameter for the conditionnal distribution of χ adjusted to match a reasonnable value of the employment rate in the current economy *)
  ];

jlim=Jmax;
(* For[j=Jmax,j≥0,j=j-1,
If[(((w[j]*(ff[j-1]-ff[j]))/(Max[ff[j],10^-10]*dw))-1-(1/Paretoskill))*(((w[Jmax]*(ff[Jmax-1]-ff[Jmax]))/(Max[ff[Jmax],10^-10]*dw))-1-(1/Paretoskill))<0,jlim=j-1;j=0];
];   (* If at the upper tail of the skill distribution, the parameter ((1-F)/(w*f)) is constrained to be a constant that equals Paretoskill, then the skill distribution is described by a Pareto distribution. The density for the upper tail then corresponds to an isolesatic function of skill, at the power -(1+(1/Paretoskill)). This command then finds the highest skill level for which the approximation of the skill distribution obtained by the kernel procedure "equals" (1+(1/Paretoskill)) i.e. for which minus the elasticity of density wrt to skill goes from above to below or from below to above 1+(1/Paretoskill)  *) *)

KPareto=ff[jlim]*((w[jlim])^(1+(1/Paretoskill)));     (* Adjust the scale parameter of the Pareto function to keep continuity *)
For[j=0,j<jlim,j=j+1,fff[j]=ff[j]];
For[j=jlim,j≤jmax,j=j+1,fff[j]=(KPareto/(((w[j])^(1+(1/Paretoskill)))));]   (* Substitute for j≥jlim, the Pareto extension of the skill distribution for its Kernel approximation  *) 
Tot=Sum[((fff[j])/(Employ[w[j]]))*dw,{j,0,jmax}];
F[-1]=0;
For[j=0,j≤jmax,j=j+1,f[j]=((fff[j])/(Tot*Employ[w[j]]));   (* unconditionnal density of skill w, indexed by j *)
F[j]=F[j-1]+f[j]*dw;    (* unconditionnal CDF of skill w, indexed by j *)
];

(* The parametrization is made for each skill level w[j] The parameters of distribution of χ conditional on w are indexed by the rank j of the skill level considered (easier for the calibration)  *)
K[χ_,j_]:= ((Exp[-α[j]+β[j]*χ])/(1+Exp[-α[j]+β[j]*χ]));       (*  CDF of χ, conditional on w, i.e. on w[j], i.e. on j *)
(*  The chosen distribution of χ is defined for χ between -∞ and χmax *)
h[u_,b_,j_]:=K[P[u,b,w[j]],j]*f[j];                                                  (* mass of employed at skill w *)

govexp=-bactual;
bcurrent=bactual;
	
For[j=0,j≤jmax,j=j+1,
hcurrent[j]=Max[(10^(-10)),h[Ucurrent[j],bactual,j]];
govexp=govexp+(TaxActual[Ycurrent[j]]+bcurrent)*hcurrent[j]*dw;
];

(*  Optimal 2nd best allocation under Maximin without extensive margin *)
ZintRawls[jmax]=0;
YintRawls[jmax+1]=2*Yint[jmax];
For[j=jmax,j≥0,j=j-1,
YintRawls[j]=((((w[j]))^(1+ϵ)))*((((((1+ϵ)/ϵ)*Max[0.01,1+((ZintRawls[j])/(hcurrent[j]))*((1+ϵ)/(ϵ*w[j]))]))^(-ϵ)));(* In case where Z is negative, the Max prevents negative values, which at power -ϵ leads to problems Hence Y is constrained to be bounded from above, at a value far above the efficient one *)
If[YintRawls[j]>YintRawls[j+1],Print["Bunching at ",j," w[j-1]=",w[j-1]," w[j]=",w[j]]];
ZintRawls[j-1]=ZintRawls[j]+hcurrent[j]*dw;
];

DVintRawls[0]=0;
For[j=0,j≤jmax-1,j=j+1,DVintRawls[j+1]=DVintRawls[j]+dw*((1+ϵ)/ϵ)*((YintRawls[j])^(((1+ϵ)/ϵ)))*((w[j])^(-((1+ϵ)/ϵ)-1))];   (* Computes the difference in C-((((Y/w)))^((1+ϵ)/ϵ)) between individuals of skill w[j] and of skill w[0]*)

bintRawls=Sum[(YintRawls[j]-((((YintRawls[j])/(w[j])))^(((1+ϵ)/ϵ)))-DVintRawls[j])*hcurrent[j]*dw,{j,0,jmax}]-govexp; (* use the government budget constraint end the optiaml condition that the least skilled among employed are  indifferent between working and not to set b *)
Print["Optimal allocation Under Maximin without extensive margin : Uw0=b=",bintRawls]

(****************************************************************************************************************************************************)
(*  Optimal 2nd best allocation under Benthamite without extensive margin *)
bint=(bintRawls/2);
λint=UC[bintRawls,0,w[0]];
bmin=0;
bmax=bintRawls;
bint=(bintRawls/2);

For[k=0,k≤kmax,k=k+1,
Adj= 1;  (*  Adj expresses the intial value of (Subscript[U, (Subscript[w, 0])]) as a fraction of its value in the laissez faire economy. A priori adj is higher than 0 *)
AdjMin=0;
AdjMax=2;
 For[kj=0,kj≤kmax,kj=kj+1, 
Zint[0]=0;
Uint[0]=Adj*bintRawls;   (* If Adj=1, the least skilled worker gets the same utility as at the Maximin optimum  *)
Zmax=0;
(* begin inner loop *)
For[j=0,j≤jmax,j=j+1,
Yint[j]=((((w[j]))^(1+ϵ)))*((((((1+ϵ)/ϵ)*Max[0.01,1+((Zint[j])/(hcurrent[j]))*((1+ϵ)/(ϵ*w[j]))*Vp[InvV[Uint[j]]]]))^(-ϵ)));(* In case where Z is negative, the Max prevents negative values, which at power -ϵ leads to problems Hence Y is constrained to be bounded from above, at a value far above the efficient one *)
Cint[j]=Γ[Uint[j],Yint[j],w[j]];
gwint[j]=((UC[Cint[j],Yint[j],w[j]])/(λint));
Uint[j+1]=Uint[j]+Uw[Cint[j],Yint[j],w[j]]*dw;Zint[j+1]=Zint[j]+(((1/(λint))-(1/(UC[Cint[j],Yint[j],w[j]])))*hcurrent[j]-((UCw[Cint[j],Yint[j],w[j]])/(UC[Cint[j],Yint[j],w[j]]))*Zint[j])*dw;
Zmax=Max[Abs[Zint[j+1]],Zmax];
];
AdjOld=Adj;
If[Zint[jmax]<0,
AdjMax=Adj;Adj=((AdjMin+Adj)/2);      ,
AdjMin=Adj;Adj=((AdjMax+Adj)/2)    ;
]
(* If Z2[jmax]<0, it means that (Subscript[Z, wmin])>0, so (Subscript[q, wmin])<0. The gvt's objective decreases by shifting up utility levels. AdjMax becomes equal to Adj and Adj is in the middle between its preceeeding value and Adjmin. Otherwise, AdjMin becomes equal to Adj and Adj is in the midle between its preceeeding value and Adjmax  *)
(* end inner loop *)
 
 ];
bold=bint;
If[Sum[(Yint[j]-Cint[j])*hcurrent[j]*dw,{j,0,jmax}]>govexp+(1-Sum[hcurrent[j]*dw,{j,0,jmax}])*bint,
bmin=bint;bint=((bint+bmax)/2),
bmax=bint;bint=((bint+bmin)/2)];  (* If the budget surplus is positive, then we re-optimize with a lower b *)
λold=λint;
λint=UC[bint,0,w[0]];
If[bmax-bmin<1,k=kmax+1];
 ];  
(*ListPlot[Table[{w[j],Zint[j]},{j,0,jmax}],Joined->True,AxesLabel->{"w","Z"}]*)
For[j=0,j<jmax,j=j+1,If[Yint[j]>Yint[j+1],Print["Bunching at ",j," w[j-1]=",w[j-1]," w[j]=",w[j]]];]

For[j=0,j≤jmax,j=j+1,
CintRawls[j]=DVintRawls[j]+bintRawls+((((YintRawls[j])/(w[j])))^(((1+ϵ)/ϵ)));];
plotRawls=ListPlot[Table[{F[j],CintRawls[j]},{j,0,jmax}],Joined->True,AxesLabel->{"Y","C"},AxesLabel->{"Y","C"},PlotStyle->{RGBColor[1,0,0]}];
plotBentham=ListPlot[Table[{F[j],Cint[j]},{j,0,jmax}],Joined->True,AxesLabel->{"Y","C"},AxesLabel->{"Y","C"},PlotStyle->{RGBColor[0,0,1]}];
Show[plotRawls,plotBentham,  PlotRange->{{0,1},{0,50000}}]

Output=Table[{F[j],CintRawls[j],Cint[j]},{j,-1,jmax}];
Export["\\\\Vboxsvr\\dd\\Économie\\Données\\ERFS_2012\\sortie_France_new3.xls",Output]

(* Affiche taux d'imposition *)
(*
For[j=0,j≤jmax,j=j+1,
Tmint[j]=(1+((UY[Cint[j],Yint[j],w[j]])/(UC[Cint[j],Yint[j],w[j]])))*100;
UintRawls[j]=DVintRawls[j]+bintRawls;
CintRawls[j]=UintRawls[j]+((((YintRawls[j])/(w[j])))^(((1+ϵ)/ϵ)));
TmintRawls[j]=(1+((UY[CintRawls[j],YintRawls[j],w[j]])/(UC[CintRawls[j],YintRawls[j],w[j]])))*100;
];
gTmint=ListPlot[Table[{Yint[j],Tmint[j]},{j,0,jmax}],Joined->True,AxesLabel->{"Y","Tm"},PlotStyle->{RGBColor[0,0,1],Dashing[{0.01,0.01}]}];
gTmintRawls=ListPlot[Table[{YintRawls[j],TmintRawls[j]},{j,0,jmax}],Joined->True,AxesLabel->{"Y","Tm"},PlotStyle->{RGBColor[1,0,0],Dashing[{0.01,0.01}]}];
gTbsurYint=ListPlot[Table[{Yint[j],Yint[j]-Cint[j]+bint},{j,0,jmax}],Joined->True,AxesLabel->{"Y","T+b"},PlotStyle->{RGBColor[0,0,1],Dashing[{0.01,0.01}]}];
gTbsurYintRawls=ListPlot[Table[{YintRawls[j],YintRawls[j]-CintRawls[j]+bintRawls},{j,0,jmax}],Joined->True,AxesLabel->{"Y","T+b"},PlotStyle->{RGBColor[1,0,0],Dashing[{0.01,0.01}]}];
Show[gTmint,gTmintRawls,PlotRange->{{0,300000},All}]
Show[gTbsurYint,gTbsurYintRawls,PlotRange->{{0,300000},{0,200000}}]
Output=Table[{F[j],CintRawls[j],Cint[j],YintRawls[j],Yint[j],TmintRawls[j], Tmint[j], YintRawls[j]-CintRawls[j]+bintRawls,Yint[j]-Cint[j]+bint},{j,-1,jmax}];
Export["\\\\Vboxsvr\\dd\\Économie\\Données\\ERFS_2012\\sortie_France_complet.xls",Output];  *)